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Abstract 

Lyapunov modes are periodic spatial perturbations of phase-space states of many-particle sys- 
tems, which are associated with the small positive or negative Lyapunov exponents. Although 
familiar for hard-particle systems in one, two, and three dimensions, they have been difficult to 
find for soft-particles. We present computer simulations for soft-disk systems in two dimensions 
and demonstrate the existence of the modes, where also Fourier-transformation methods are em- 
ployed. We discuss some of their properties in comparison with equivalent hard-disk results. The 
whole range of densities corresponding to fluids is considered. We show that it is not possible to 
represent the modes by a two-dimensional vector field of the position perturbations alone (as is the 
case for hard disks), but that the momentum perturbations are simultaneously required for their 
characterization. 

PACS numbers: 05.45.PqNumerical simulation of chaotic systems, 05.20.-yClassical statistical mechanics, 
47.35.+iHydrodynamic waves 
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I. INTRODUCTION 



For the last 50 years molecular dynamics simulations have decisively contributed to our 
understanding of the structure and dynamics of simple fluids and solids [1]. More recently, 
also the concepts of dynamical systems theory have been applied to study the tangent- 
space dynamics of such systems [2j. Of particular interest is the extreme sensitivity of 
the phase-space evolution to small perturbations. On average, such perturbations grow, or 
shrink, exponentially with time, which may be characterized by a set of rate constants, the 
Lyapunov exponents. The whole set of exponents is referred to as the Lyapunov spectrum. 
This instability is at the heart of the ergodic and mixing properties of a fluid and offers a 
new tool for the study of the microscopic dynamics. In particular, it was recognized very 
early that there is a close connection with the classical transport properties of systems in 
nonequilibrium stationary states Q]. For fluids in thermodynamic equilibrium, an analysis 
of the Lyapunov instability is expected to provide an unbiased expansion of the dynamics 
into events, which, in favorable cases, may be associated with qualitatively different degrees 

inear molecules 14, or with the intra- 



of freedom, such as the translation and rotation of linear molecules 
molecular rotation around specific chemical bonds 



|5|. 



Since the pioneering work of Bernal |(| with steel balls, hard disks have been considered 
the simplest model for a "real" fluid. With respect to the structure, they serve as a reference 
system for highly-successful perturbation theories of liquids b l |gl_ - Recently we studied the 
Lyapunov instability of such a model and found 



1. that the slowly-growing and decaying perturbations associated with the non- vanishing 
Lyapunov exponents closest to zero may be represented as periodic vector fields co- 
herently spread out over the physical space and with well-defined wave vectors k. 
Because of their similarity with the classical modes of fluctuating hydrodynamics we 
refer to them as Lyapunov modes. Depending on the boundary conditions, the respec- 
tive exponents are degenerate, and the spectrum has a step-like appearance in that 
regime. 

2. that the fast-growing or decaying perturbations are localized in space, and the number 
of particles actively contributing at any instant of time vanishes in the thermodynamic 
limit. 
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Experimentally, the Lyapunov modes were found for hard-particle systems in one, two 
and three dimensions and for various boundary conditions, provided that the system size, L, 
in at least one direction is large enough for the discrete particles to generate a recognizable 
wave-like pattern with a wave number k. Their theoretical understanding is based on the 
spontaneous breaking of the translational symmetry of the zero modes (k = 0), which are 
associated with the vanishing Lyapunov exponents and are a consequence of the conservation 



laws in the system |l3, 
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Based on this experimental and theoretical evidence for hard-particle systems, it is natural 
to expect that the Lyapunov modes are a general feature of many-body systems with short- 
range interactions and that the details of the pair potential should not matter. However, 
already our very first simulations with soft disks in two dimensions 0, H revealed 
a much more complicated scenario. Whereas the spatial localization of the fast-growing or 
decaying perturbations could be easily verified, the mode structure for the slowly-growing or 
decaying perturbations was elusive and could, at first, not be unambiguously detected. Here 
we demonstrate that the modes for soft sphere fluids do indeed exist. However, sophisticated 
Fourier-transformation techniques are needed to prove their existence. Interestingly, the 
degeneracy of the Lyapunov exponents and, hence, the step structure of the spectrum so 
familiar from the hard-disk case is recovered, but only for low densities. This suggests that 
kinetic theory is a proper theoretical framework in that case. For intermediate and large- 
density soft-disk fluids the degeneracy disappears. We do not have an explanation for this 
fact. Most recently, Lyapunov modes were demonstrated by Radons and Yang 0, |2l| for 
one-dimensional Lennard- Jones fluids at low temperatures and densities. 

In this work we analyze the Lyapunov instability of repulsive Weeks- Chandler- Anderson 
(WCA) disks in two dimensions and compare it to analogous hard-disk results. The main 
difficulty is the large box size required for the modes to develop and, hence, the large number 
of particles, N. It requires parallel programming for the computation of the dynamics both in 
phase and tangent space and for the re-normalization of the perturbation vectors according 
to the classical algorithms of Benettin et al. and Shimada et al. |23j. In Sec. ITT1 we 
characterize the systems and summarize the numerical methods used. In Sec. IHII we discuss 
the surprising differences found between the Lyapunov spectra for the soft- and hard-particle 
fluids, particularly at intermediate and large densities. Sec JIVI is devoted to an analysis of 
the zero modes associated with the vanishing Lyapunov exponents. The Lyapunov modes 
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are analyzed in Sees. IVland lVll We close with some concluding remarks in Sec. IVHI 



II. CHARACTERIZATION OF THE SYSTEM AND METHODOLOGY 

We consider N disks with equal mass m in a two-dimensional rectangular box with 
extensions L x and L y and aspect ration A = L y /L x . Periodic boundary conditions are used 
throughout. The particles interact with a smooth repulsive potential, for which we consider 
two cases: 

i) a Weeks-Chandler-Anderson potential, 

12 



>WCA 




r > 2 1 / 6 ct 



with a force cutoff at r = 2 1 / 6 cr. As usual, a and e are interaction-range and energy 
parameters, respectively, which are ultimately set to unity. Actually, such a potential is not 
the best choice, since the forces and their derivatives (required for the dynamics in tangent- 
space) are not continuous at the cutoff. This introduces additional noise into the simulation 
and violates the conservation laws. It also affects the computation of the Lyapunov spectrum. 
To avoid this problem we sometimes also use 
ii) a power-law potential, 



lOOe 
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r < a 

(2) 

0, r > a 



which looks very similar but does not suffer from this deficiency. We mention, however, that 
the results discussed below turn out to be insensitive to this weakness of the WCA potential. 
For reasons of comparison we also consider 
iii) a hard-disk potential, 

{oo, r < a 

(3) 
0, r > a 

The phase space, X, has AN dimensions, and a phase point r £ X is given by the 4iV- 
dimensional vector F(t) = {qi,pi]i = 1,...,N}, where $ and pi, denote the respective 
positions and momenta of the disks. It evolves according to the time-reversible motion 
equations 

f = F(T), (4) 



which are conveniently written as a system of first order differential equations. Here, F(T) = 
{ m' — f^' * — !)•••) ^} follows from Hamilton's equations, where $ = £\ ^j>i 0(1?? — 
is the total potential energy. 

Any infinitesimal perturbation 5Ti = {8qf , 8pf , i = 1, . . . , N} lies in the tangent space 
TX, tangent to the manifold X at the phase point T(t). Here, 8q\ , and 8pp are two- 
dimensional vectors and denote the position and momentum perturbations contributed by 
particle i. STi evolves according to the linearized equations of motion 

dF 

5T l = ^ r - 5T l . (5) 
dTi K J 

Oseledec's theorem |24j assures us that there exist 4N orthonormal initial tangent vectors 
5Ti(0),l = 1, . . . ,4iV, whose norm grows or shrinks exponentially with time such that the 
Lyapunov exponents 

A -=fcT ln i^' '=i.-.*". («) 



exist. However, in the course of time these vectors would all get exponentially close to 
the most-unstable direction and diverge, since the unconstrained flow in tangent space does 
not preserve the orthogonality of these vectors. In the classical algorithms of Benettin et 



al. |22|, and Shimada et al. |23[ also used in this work, this difficulty is circumvented by 
replacing these vectors by a set of modified tangent vectors, 8i,l = 1, . . . ,4iV, which are 
periodically re-orthonormalized with a Gram-Schmidt procedure. The Lyapunov exponents 
are determined from the time-averaged renormalization factors. The vectors 8i represent the 
perturbations associated with A; and are the objects analyzed below. Henceforth we use the 
notation 8qf \ and 8pf' also for the individual particle contributions to the ortho-normalized 
vectors, 8i = {8q\ , 8pf , 2 = 1,..., iV}. For later use, we introduce also the squared norm 
in the single-particle \i space, 

2 / , rt \ 2 



7<" 



(««?>) + {Spf) . (7) 



This quantity is bounded, < 7^ < 1, and obeys the sum rule ^^l 1 7f' > = 1 for any I. 
Since the equations (0) are linear in 8Ti, the quantities 7^ indicate, how the activity for 
perturbation growth, as measured by A/, is distributed over the particles at any instant of 
time. 
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For the hard-disk systems the motion equations (@J) and (J^J) need to be generalized to 
include also the collision maps due to the instantaneous particle collisions. For algorithmic 



details we refer to our previous work [9 



25, 



The exponents are taken to be ordered by size, A^ > Aj+i, where the index I numbers 
the exponents. According to the conjugate pairing rule for symplectic systems 0Q> the 
exponents appear in pairs such that the respective pair sums vanish, A; + Xm+i-i = 0. Only 
the positive half of the spectrum needs to be calculated. Furthermore, six of the exponents, 
{A 2 AT-2<z<2A r +3} ) vanish as a consequence of the constraints imposed by the conservation 
of energy, momentum, center of mass (for the tangent-space dynamics only), and of the 
non-exponential time evolution of the perturbation vector in the direction of the flow. 

In the soft-disk case we use reduced units for which the disk diameter a, the particles' 
mass m, the potential parameter e, and Boltzmann's constant fee are all set to unity. As 
usual, the temperature is defined according to (K) = (^z^ = (N — l)k B T, where K is 
the total kinetic energy and (...) is a time average. For hard disks if is a constant of the 
motion, and K/N is taken as the unit of energy. Since Lyapunov exponents for hard disks 
strictly scale with the particle velocities and, hence, the square root of the temperature, all 
comparisons between soft and hard disks below are for equal temperatures, T — 1. 

The computation of full Lyapunov spectra for soft particles is a numerical challenge, since 
the number of simultaneously integrated differential equations increases with the square of 
the particle number. Therefore, parallel processing with up to 6 nodes is used, both for the 
integration and for the Gram-Schmidt re-orthonormalization procedure j^. 



III. LYAPUNOV SPECTRA FOR SOFT-DISK FLUIDS 



A. Phenomenology 

We start our comparison of the spectra for soft (WCA) and hard-disks with low-density 
gases, p = N/L x L y = 0.2, at a common temperature T — 1. The system consists of N = 80 
particles in a rather elongated rectangular simulation box, L x = 100, L y = 4; A = L y /L x = 
0.04, with periodic boundary conditions. An inspection of the upper panel of Fig. shows, 
as expected, that the two systems have similar spectra and, thus, similar chaotic properties: 
the maximum exponents, Ai, and the shapes of the spectra are reasonably close. For even 
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lower densities (not shown) the agreement is even better. Most important, however, is the 
observation that the step structure due to the degeneracy of the small Lyap unov exponents 



is observed both for the hard disks, for which it is well known to exist [13J, |25j, |30j, |3l|, and 
for the soft-disk gases at low densities. We conclude from this comparison in Fig. that 
the Lyapunov modes exist for low-density soft disks. In view of the smallness of L y and 
the quasi one- dimensional nature of the system, these modes may have only wave vectors 
parallel to the x axis of the box. 

This general picture also persists for intermediate gas densities, p = 0.4, as shown in the 
lower panel of Fig. ^ The box size is the same as before, but the number of particles (and, as 
a consequence, also the number of exponents) is doubled. The spectral shapes for the WCA 
and hard-disk fluids, and the maximum exponents Ai in particular, become significantly 
different. Chaos, measured by Ai, is obviously enhanced for the WCA particles with a finite 
collision time as compared to the instantaneously colliding hard disks. Interestingly, the 
Kolmogorov-Sinai (or dynamical) entropy hxs, which is equal to the sum of all positive 
exponents differs less due to the compensating effect of the intermediate exponents 

in the range 0.2 < 1/2N < 0.8. This may be verified in Fig. HI below. 

The systems of Fig. ^ are quasi one-dimensional and do not represent bulk fluids. We 
show in Fig. El analogous spectra for bulk systems in a rectangular simulation box with 
aspect ratio A = 0.6, and containing A" = 375 particles. The density varies between 0.1 and 
0.4 and is specified by the labels. For all spectra the lowest step is clearly discernible, but 
is progressively less pronounced if the density is increased. For p = 0.1, the four smallest 
exponents (744 < I < 747) are identified as belonging to longitudinal modes (see below), 
the two next-larger exponents (742 < I < 743) to transverse modes. In Figure H3 we compare 
a WCA system to a hard-disk system in a square simulation box at a density p = 0.4 and 
A" = 400 particles. The length of the simulation box in every direction is L x = L y = 31.62. 
The step structure that is so prominent in the hard disk spectrum vanishes altogether for the 
soft potential case. If the density of such a square system with A" = 400 particles is reduced 
by increasing the size of the simulation box, the step structure with the same degeneracy 
reappears. 
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FIG. 1: Lyapunov spectra for soft (WCA) and hard (HD) disks at low (p = 0.2, upper panel) and 
intermediate (p = 0.4, lower panel) densities. The very elongated rectangular simulation box is the 
same in both cases (L x = 100, L y = 4, A = 0.04), such that N = 80 (top) and N = 160 (bottom) 
disks are involved, respectively. The temperature T = 1. The Lyapunov index I is normalized by 
2N. Only the positive branches of the spectra are shown. Although the spectra are only defined 
for integer values of I, smooth lines are drawn for clarity. In the inset a magnified view of the 
small-exponents regime is shown, where I is not normalized. 

B. Measures for chaos 



The maximum Lyapunov exponent, Ai, and the Kolmogorov- Sinai entropy, hxs, are 
generally accepted as measures for chaos. The latter corresponds to the rate of information 
generated by the dynamics and, according to Pesin's theorem 32, is equal to the sum 
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I/2N 

FIG. 2: Lyapunov spectra for WCA-disk fluids with N = 375 particles in a periodic box with fixed 
aspect ratio A = 0.6. Only the small exponents are shown as a function of the reduced index 1/2N. 
The density is indicated by the labels. 
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FIG. 3: Lyapunov spectrum for WCA and hard-disk fluids with N = 400 particles, a density 
p = 0.4, and L x = L y = 31.62. The temperature T = 1. The Lyapunov index I is normalized by 
2N. The inset shows the lower part of the spectra where the step structure is prominent for the 
hard disks but totally absent for the WCA particles. For this part, the index is not normalized. 

of the positive Lyapunov exponents, hxs = J2\ t >o ^43- IHwe compare the isothermal 

density dependences of these quantities for 375-disk WCA systems with corresponding hard- 
disk results. Also included in this set of figures is a comparison for the smallest positive 
exponent, A27V-3, which is directly affected by the mechanism generating the Lyapunov mode 
with the longest-possible wave length (if it exists). 
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FIG. 4: Isothermal density dependence of the maximum Lyapunov exponent, Ai (top), of the 
smallest positive exponent, \2N-3 (middle), and of the Kolmogorov-Sinai entropy per particle, 
hics/N (bottom), for hard and soft-disk systems as a function of the density. The particle number, 
N = 375, the aspect ratio, A = 0.6, of the periodic box, and the temperature, T = 1 are held fixed. 
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For very dilute hard-disk gases, Ai and hxs/N have been estimated from kinetic theory by 
Dorfman, van Beijeren and van Zon j^j], who have provided expansions to leading orders in 
p, Ax = Ap(— In p + B) + . . . , and hxs/N = Ap(— In p + B) + . . . , with explicit expressions 
for the constants A, B and A, B. These predictions were very successfully confirmed by 
computer simulations j3, Computer simulations by Dellago, Posch and Hoover also 
provided hard-disk results over the whole range of densities j^j]. Ai and hxs/N were found to 
increase monotonously with p with the exception of a loop in the density range between the 
freezing point of the fluid ( pf = 0.88 js?!) and the melting point of the solid (p m = 0.91 ^\). 
These loops disappear if Ai and hxs/N, are plotted as a function of the collision frequency 
v instead of the density j25|. If p approaches the close-packed density p — 1.1547/<t 2 , both 
Ai and h KS /N diverge as a consequence of the divergence of v. 

The results for the WCA disks reported here are for 0.1 < p < 0.8, which covers almost 
he whole fluid range. Simulations for rarefied gases and for solids are currently under way 
3f^ . Not unexpectedly, we infer from the top panel of Fig. HJthat Ai-WCA approaches Ai-HD 
for small densities p < 0.1. However, The density dependence for larger p is qualitatively 
different. Ai approaches a maximum near the density corresponding to the fluid-to-solid 
phase transition, which confirms our previous results for Lennard- Jones fluids and solids 
39]. The collective dynamics at such a transition causes maximum chaos in phase space. 
In the solid regime (not included in Fig. 0]) Ai drops with increasing density, which is easily 
understandable. 

The smallest positive exponent A2at-3 for soft-disk gases follows the hard-disk density 
behavior more closely than Ai for lower densities as inferred from the middle panel of Fig. 
0] This means that the associated collective dynamics is affected less by the details of the 
pair potential. Significant deviations start to occur for p > 0.3. For }iks/N the agreement 
between soft and hard disks extends even to a larger range of densities, as may be seen in 
the bottom panel of Fig. HJ However, this observation is deceptive and does not necessarily 
signal a similar dynamics. It is a consequence of significant cancellation due to the exponents 
intermediate between the maximum and the small exponents of the spectra such as in the 
bottom panel of Fig. [TJ 
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C. Localization of tangent-space perturbations in physical space 



One may interpret the maximum (minimum) Lyapunov exponent as the rate constant 
for the fastest growth (decay) of a phase-space perturbation. Thus, it is dominated by the 
fastest dynamical events, binary collisions in the case of particle fluids. Therefore, it does not 
come as a surprise that the associated tangent vector has components which are strongly 
localized in physical space |4j. Similar observations for other spatially-extended systems 
have been made before by various authors 0, 41, 12, Q- With the help of a particular 



measure for the localization 



12|, we could show that for both hard and soft disk systems 



the localization persists even in the thermodynamic limit, such that the fraction of tangent- 
vector components contributing to the generation of Ai at any instant of time converges to 
zero with N — > oo The localization becomes gradually worse for larger Lyapunov 

indices I > 1, until it ceases to exist and (almost) all particles collectively contribute to the 
perturbations associated with the smallest Lyapunov exponents, for which coherent modes 
are known (or believed) to exist Isi). 

Here, we adopt another entropy-based localization measure due to Taniguchi and Mor- 
riss [44] which was successfully applied to quasi-onedimensional hard-disk gases. Recalling 
the definition of the individual particle contributions jf' to Sf(= 1) in Eq. (J7J), one may 
introduce an entropy-like quantity 

N 

^ ) = -E^ (*)W°(*)>, (8) 

1=1 

where (...) denotes a time average. The number 

W w = exp (S®) (9) 

may be taken as a measure for the spatial localization of the perturbation 5i associated with 
Az, 1 < I < 2N. It is bounded according to 1 < W"' < N. The lower bound, 1, indicates 
the most-localized state with only one particle contributing, and the upper bound N signals 
uniform contributions by all of the particles {7P = l/N,i = 1, . . . ,N}. We shall refer to 
the set of normalized localization parameters {W^/NJ = 1, . . . ,2N} as the localization 
spectrum. 

In Fig. El we compare the localization spectra for the same hard-disk and soft-disk sys- 
tems for which the Lyapunov spectra are given in Fig. ^ (L x = 100, L y = 4, A = 0.04). 
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FIG. 5: Dependence of the normalized localization parameters, W^'/N, on the normalized Lya- 
punov index 1/2N. The corresponding Lyapunov spectra are shown in Fig. ^ The soft and hard 
disks are indicated by the labels and distinguished by color. The system is a rectangular box 
(L x = 100, L y = 4) with periodic boundaries. Top: p = 0.2, N = 80. Bottom: p = 0.4, N = 160. 
The temperature T = 1. The insets provide a magnified view of the mode regime. 

For the lower-density gases in the top panel (p = 0.2, iV = 80), the localization spectra for 
soft and hard disks are almost indistinguishable. They indicate strong localization for the 
perturbations belonging to the large exponents (small /), and collective behavior for large /. 
Of particular interest is the comb-like structure also magnified in the inset, which is a con- 
sequence of the Lyapunov modes (which will be discussed in more detail below). Stationary 
transverse modes are hardly affected by the time averaging involved in the computation of 
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the localization measure and lead to large values of W^/N, indicating strong collectivity. 
The propagating LP-modes ^ , however, are characterized by a significantly- reduced value 
of W^/N. For intermediate gas densities in the bottom panel (p = 0.4, N = 160), differ- 
ences between the hard and soft-disk results become apparent. Most prominently, the comb 
structure in the localization spectrum of the WCA system has mostly disappeared, although 
the Lyapunov spectrum clearly displays steps in the small-exponent regime, which is a clear 
indication for the existence of modes. 



IV. ZERO MODES 

The dynamics in phase space and tangent space is strongly affected by the inherent 
symmetries of a system, i. e. infinitesimal transformations leaving the equations of motion 
invariant. They are intimately connected with the conservation laws obeyed by the dynamics. 
If these transformations act as infinitesimal perturbations of the initial conditions, the latter 
do not grow/shrink exponentially (but at most linearly) with time and give rise to as many 
vanishing Lyapunov exponents. The generators of these transformations are taken as unit 
vectors in tangent space, which point into the direction of the respective perturbation and are 
called zero modes. They span an invariant subspace AA(T) of the tangent space at any phase 
point T, which is referred to as the zero subspace. They are essential for an understanding 
of the Lyapunov modes dealt with in the following. 

For a planar system of hard or soft disks, there are six symmetry-related perturbations 
leading to non-exponential growth lJ 3| and, hence, six vanishing Lyapunov exponents: 



1. Homogeneous translation in x and y directions, with generators e\ and e2 explicitly 
given below. It is a consequence of center-of-mass conservation, which is actually not 
obeyed for systems with periodic boundary conditions, but, nevertheless, still holds 
for the linearized motion in tangent space, which does not recognize periodic box 
boundaries for the dynamics of the perturbations. 

2. Homogeneous momentum perturbation in x and y directions due to Galilei invari- 
ance, with generators and given below. It may be viewed as a consequence of 
momentum conservation. 

3. Simultaneous momentum and force rescaling according to a generator e§ as given 
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below. It is a consequence of energy conservation and, ultimately, of the homogeneity 
of time. 



4. Homogeneous time shift with generator corresponding to non-exponential growth 
in the direction of the phase flow and also a consequence of the homogeneity of time. 

These generators follow from the constraints the conservation equations impose on the 
dynamics in phase space. The latter constitute hyper-surfaces, Y2iLi x i = cons t, YliLiyi = 
const for the center of mass (see the comment above), YliLxPx,i = 0? "Y^H=iPy,i = f° r 
the conserved momentum, and Y2iLi(Pxi ^~Pyi)/^ + ^ = F for the conserved energy. The 
gradients perpendicular to these hyper- surfaces provide directions of non-exponential phase- 
space growth and, hence, the generating vectors e%, . . . , e$ If we use the explicit notation 



ST = {Sxi, Syi, 5p Xji , 8py }i ; i = 1, . . . , N} 
for an arbitrary tangent vector, one obtains 

e x = (l/iV){l,0,0,0;z = l,...,iV} (10) 

e 2 = (l/iV){0,l,0,0;z = l,...,iV} (11) 

e 3 = (l/JV){0,0,l,0;i = l,...,iV} (12) 

e 4 = (l/JV){0,0,0,l;z = l,...,iV} (13) 

^5 ^{ F X) ij ^y,ii Px,ii Py,ii ^ !>•••) ^} (-^) 

e 6 = a{p Xti ,p y> i,]F x> i,F y)i ;i = l,...,N}, (15) 



where a is a normalizing factor. Each vector has AN dimensions. e\ to e% are orthonormal 
and form a natural basis for the invariant zero space. In the following we consider the 
expansion of the six orthonormal tangent vectors S 2 n-2, ■ ■ ■ ,^2N+3, responsible for the six 
vanishing exponents in the simulation, in that basis, 

6 

<W-3+i = ^2 ai 'i e i ' 1 = • ' ' ' 6 ' ( 16 ) 

with projection cosines, ay = (5i • ej). Since all vectors involved are of unit length, ay may 
either be interpreted as the projection of Son-z+i onto ej, or vice versa. 

For hard disk fluids one can easily show |13J| that the tangent vectors 82N-2, ■ ■ ■ > 82N+3 are 
fully contained in Af(T). The projection cosines strictly obey the sum rules 

6 6 

^(a M ) 2 = l , ^(a M ) 2 = l, 
i=i 1=1 
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TABLE I: Projection cosines, aij, according to Eq. (|16|) for N = A soft particles in a periodic box 
and interacting with the power-law potential (J5J). The density is low, p = 0.1, to isolate individual 
collision events. Two instants are considered: 

A: At time to just before a collision of particles 1 and 2: 



I 


a/,i a ij2 «/,3 a^ 4 


a/,5 


a/,6 


£, 6 =i Hi) 2 


1 


-1.00000 0.00000 0.00000 0.00000 


0.00000 


0.00000 


1.00000 


2 


0.00000 0.21120 0.00000 0.00000 


-0.00001 


-0.97734 


0.99980 


3 


0.00000 0.97744 0.00000 0.00000 


0.00000 


0.21118 


0.99999 


4 


0.00000 0.00000 0.95954 0.00000 


0.28155 


0.00000 


0.99998 


5 


0.00000 0.00000 0.28158 0.00000 


-0.95944 


0.00001 


0.99981 


6 


0.00000 0.00000 0.00000 1.00000 


0.00000 


0.00000 


1.00000 


Eli K;) 2 


1.00000 1.00000 1.00000 1.00000 


0.99979 


0.99979 





B: At time to + 0.11 of closest approach of particles 1 and 2: 



I 


a/,i aj i2 a/ j3 az j4 


a/,5 


a/,6 


£?=i (ay) 2 


1 


-1.00000 0.00000 0.00000 0.00000 


0.00000 


0.00000 


1.00000 


2 


0.00000 0.14132 0.00000 0.00000 


0.00000 


-0.08601 


0.02737 


3 


0.00000 0.98996 0.00000 0.00000 


0.00000 


0.01228 


0.98018 


4 


0.00000 0.00000 0.98172 0.00000 


0.01654 


0.00000 


0.96405 


5 


0.00000 0.00000 0.19033 0.00000 


-0.08529 


0.00000 


0.04350 


6 


0.00000 0.00000 0.00000 1.00000 


0.00000 


0.00000 


1.00000 


Eli K;) 2 


1.00000 1.00000 1.00000 1.00000 


0.00755 


0.00755 





at all times, if the total momentum vanishes and if K — N, as is always the case in our 
simulations. Thus, Span(5 2 /v-2, • • • , S2N+3) = Span(ei, . . . , e^). Note that there are no forces 
acting on the particles in this case, and the particles are moving on straight lines between 
successive instantaneous collisions. 

For soft disk systems, however, the situation is more complicated. In Table |l] we list 
instantaneous matrix elements aij for a simple example, N = 4 particles interacting with 
the particularly-smooth repulsive potential of Eq. El The density, p = 0.1, is low enough 
such that isolated binary collisions may be easily identified. Two instances are considered: 
one just at the beginning of a collision of particles 1 and 2 (upper part of the table), and a 
time half way through this collision (lower part of the table). Considering the columns first, 
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D, D_ 




FIG. 6: Caricature of the high-dimensional tangent space for N particles. The orthogonal zero 
modes e% and &2 always point into the subspace T> + = Span(52jV-2 5 $2N-i, ^2n)- Similarly, the zero 
modes and e 4 always point into the subspace D_ = Span(J2Af+i, <^2Af+2, ^2^+3)- The subspaces 
V + and D_ are represented by the thick arrows. The projection angle of the zero mode 65 into T> + 
agrees with that of e§ into P_ and generally differs from zero whenever a collision occurs anywhere 
in the system. 

the squared sums Ym=i ( a i,j) 2 always add to unity for j — 1, . . . , 4, which means that e\ to e 4 
are fully contained in the subspace Span(<5 2A r_ 2 , • • • , (W+3) of the tangent space. The same 
is not true, however, for and e^, which contain in their definition the instantaneous forces 
acting on the colliding particles. If no collisions take place anywhere in the system (upper 
part), the subspaces Span(<5 2 7v-2, • • • , 82N+3) and Span(ei, . . . , e^) are nearly the same, but 
not quite, as the squared projection-cosine sums indicate. The moment a collision occurs, 
the sums Y^n=i ( a i,j) 2 f° r J = 5 and 6 may almost vanish, as happens in the bottom part 
of the table, and the vectors e 5 and e 6 become nearly orthogonal to Span(5 2 Ar_2, • • • , <5 2 tv + 3). 
In Figure |H1 an attempt is made to illustrate this relationship for such a high-dimensional 
tangent space. For systems containing many particles and/or for larger densities, there 
will always be at least one collision in progress, and the horizontal sums, Y^j=i ( a i,j) 2 f° r 
I — 1, ... ,6, and the vertical sums, Ym=i i a i,i) 2 f° r J = 5 and 6, fluctuate and assume values 
significantly smaller than unity. The subspaces Span(5 2 Ar-2, • • • , = © 22- (see the 

definition in Fig. EJ) and Af = Span(ei, . . . , e 6 ) do not agree in general. 
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V. LYAPUNOV MODES 



Lyapunov modes are periodic spatial patterns observed for the perturbations associated 
with the small positive and negative Lyapunov exponents. They are characterized by wave 
vectors k with wave number 




k n x ,n y = \ (-j-n x ) +[-j-n y ) ; n x ,n y = 0,1, .. . , (17) 



where a rectangular box with periodic boundaries is assumed, and where n x and n y denote 
the number of nodes parallel to x and v, respectively. For modifications due to other 
boundary conditions we refer to Refs. [13|, |45( . Experimentally, Lyapunov modes were 
observed for hard particle systems in one, two, and three dimensions |9|, |ll|, [12 



or hard planar dumbbells [Q, |12j, |46| , and, most recently, for one-dimensional soft particles 
20l |21| . Theoretically, they are interpreted as periodic modulations (k ^ 0) of the zero 
modes, to which they converge for k — > 0. Since a modulation, for k > 0, involves the 
breaking of a continue symmetry teanslationa. symmetry of the zero modes), they have 
been identified as Goldstone modes |12lj, analogous to the familiar hydrodynamic modes and 
phonons. For the computation of the associated wave- vector dependent Lyapunov exponents 
governing the exponential time evolution of the Lyapunov modes, random matrices jl^L I30I ] , 
periodic-orbit expansion j^|, and, most successfully, kinetic theory have been employed 

HHQ. 

For hard-disk systems the representation of the modes is simplified by the fact that, for 
large N, the position perturbations, {5qf\i = 1,...,N}, and momentum perturbations, 
{6pf\ i = 1, . . . , N}, may be viewed as two-dimensional vector fields (fiq\x, y) and (fp\x, y), 
respectively, which turn out to be nearly parallel (for A > 0), or anti-parallel (for A < 0) 
1 1 31 ] . To illustrate this point, we plot in Fig. [7| 



<oos(e,)> = ( ^/r^r;/^ > • w 



where B; is the angle between the 2iV-dimensional vectors comprising the position and 
momentum perturbations of 5i, when they are viewed in the same 2iV-dimensional space. 
As always, (...) denotes a time average. Only the indices / < 2N — 2 corresponding to 
positive exponents are considered. G; (upper blue curve) nearly vanishes in the mode regime 
(close to the right-hand boundary of that figure). For the negative exponents (not shown) 
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FIG. 7: Time-averaged values for cos(0/), as defined in Eq. (|18|) . for the indices 1 < I < 2N — 3 
corresponding to the positive Lyapunov exponents of a spectrum. The systems consist of N = 375 
particles at a density p = 0.4 in a rectangular periodic box with aspect ratio A = 0.6; Upper blue 
points: hard disks (HD); Lower red points: soft WCA disks. In both cases, the temperature is 
unity. 

the angle approaches it. For large-enough A" and small positive A/, the individual particle 
contributions behave as 



where C\ is a positive number, which is almost the same for all particles i. Once the 5q\ 
are known, the 8pf may be obtained from this relation. For a characterization of the modes, 
it suffices to consider only the position perturbations, which are interpreted as a vector field, 
tpq\x,y), over the simulation cell jl^ ]. 

For soft-spheres the situation is more complicated. For low-density systems, cos(0/) 
behaves as for hard disks, which is expected in view of the similarity of the Lyapunov 
spectra. For dense fluids, however, this quantity behaves qualitatively different and does 
not converge to unity but to zero in the interesting mode regime. This is demonstrated in 
Fig. [7| by the lower red curve, which is for A" = 375 soft WCA disks at the same density 
p = 0.4. The larger the collision frequency, the larger the deviations of (cos(G/)) from 1 
(for A/ > 0), or -1 (for A/ < 0) may become. This result indicates that a representation of a 
Lyapunov mode in terms of a single (periodic) vector field of the position (or momentum) 
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perturbations of the particles is possible only for low-density systems. For larger densities 
and/or larger systems with many collisional events taking place at the same time, one needs 
to simultaneously consider all the perturbation components 5xi, 6z/i, Sp Xti , 8p yti of a particle. 

To understand this unexpected behavior in more detail, we display in the top panel of Fig. 
[SJ the time dependences of cos(O^) belonging to the maximum (/ = 1) and to the smallest- 
positive (I = 2N — 3 = 5) exponents of a four-particle system in a square periodic box with 
a density p = 0.1. During a streaming phase without collisions, the 2iV-dimensional vectors 
8q® and 8p^ tend to become parallel as required by the linearized free-flight equations 
in tangent space. However, this process is disrupted by a collision, which reduces (and 
in some cases even reverses the sign of) cos(@;). In the following streaming phase, i.e. 
forward in time, Qi(t) relaxes towards zero. This suggests that the numerical time evolution 
of cos(Oi) and, hence, of 5i is not invariant with respect to time reversal. This is indeed 
the case as is demonstrated in the lower panel of Fig. |H1 To construct this figure, the 
phase space trajectory was stored for another 10000 time units and, after a time-reversal 
transformation, {qi — > q^p% — > —p%,i = 1,...,N}, was consecutively used - in reversed 
order - as the reference trajectory for the reversed tangent-space dynamics. In the lower 
panel we show cos(Gig) and cos(Oi2), which belong to the minimum and largest-negative 
exponents, of the time-reversed dynamics, respectively, and which should be compared to 
cos(Gi) and cos(G>5) for the forward evolution. Although the same collisions are involved, the 
curves look totally different. This confirms previous results concerning the lack of symmetry 
for the forward and backward time evolution of the Gram-Schmidt orthonormalized tangent 
vectors Si for Lyapunov-unstable systems Furthermore, we have numerically verified 

that cos(6/(t)) = cos(©4jv-z+i(£)) is always obeyed, both forward and backward in time. 

VI. FOURIER-TRANSFORM ANALYSIS 

Because of numerical fluctuations, a simple inspection of the vector field (p q (x, y) very 
often is not sufficient to unambiguously determine the wave vector of a mode, particularly 
for the smaller than the maximum-possible wave lengths. Therefore, Fourier transformation 
methods have been used where we regard ip a (x,y) as a spatial distribution 



N 




(19) 



i=i 
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FIG. 8: Time dependence of cos(@;) both in the time-forward (upper panel) and time-backward 
directions (lower panel) for four WCA particles in a square periodic box. The density p = 0.1. The 
respective indices, I, are indicated by the labels. For comparison, an analogous hard-disk result for 
I = 1 is also shown by the blue curve in the upper panel. 

The Oj are identified, for example, with 5xf' or with 



7i 



' : ] ~ f(5xf^) 2 + (5yf^) 2 + (5p x 1 ^) 2 + (Spyl) 2 ) ■ In view of the periodicity of the box, 



the Fourier coefficients have wave numbers given by Eq. p7jl . They are computed from 

' r d 2 qe k -iXa(q) = - r ^-Y. a * ek " }l - 



Xa{k) 



L X Ly 



L X Ly 



(20) 



The power spectrum is defined by 



Pa{k) = Xa(k) X a(-k). (21) 

For the a; identified above, the power spectra are denoted by P x {k) and P^^ik), respectively. 
We have also applied the algorithm for unequally-spaced points by Lomb , , suitably 
generalized to two-dimensional transforms. 
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FIG. 9: Power spectra P qp {k) for the transverse modes T(1,0) to T(6,0) with I = 317, 311, 305, 
299, 293, and 287, for N = 160 WCA particles in an elongated periodic box, L x = 100, L y = 4. 
The density p = 0.4. The Lyapunov spectrum for this system is shown in the lower panel of Fig. ^ 

As a first example we show in FigureElthe power spectra Py/2(fc) for successive transverse 
modes, indexed by / = 317, 311, 305, 299, 293, and 287, for the WCA-system described in 
the lower panel of Fig. ^ They correspond to the successive wave numbers ki$ to k^o 
as predicted by Eq. (|17|). In a second and less-trivial example, we show in Fig. ITU1 the 
Lyapunov spectra (left panel) of WCA particles in a fixed elongated box, L x = 100, L y = 4, 
for various densities varying from p — 0.1 to p — 0.7 as indicated by the labels. In the 
right panel the power spectra P_,i/2(&) for the modes corresponding to the smallest positive 
exponent, A; = 2at_3, of each spectrum are shown. Since L x is the same in all cases, all power 
spectra have a peak at the allowed wave number ki t o = 0.063. This peak is also well resolved 
for large densities, for which the step structure in the spectrum is blurred. 

VII. CONCLUDING REMARKS 

In the foregoing sections we have demonstrated the existence of Lyapunov modes in 
soft-disk fluids with the help of various indicators such as the localization measure in Sec. 
1111 CI and the Fourier analysis in the previous section. However, the classification and 
characterization of the modes is more complicated than for hard disks. This is demonstrated 
in Fig. ^2 where the position perturbations 5x[ m (bottom) and SyJ 46 (top) are plotted 
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FIG. 10: WCA particles in a fixed periodic box, L x = 100, L y = 4, for various densities as indicated 
by the labels. The temperature is unity. Upper panel: Lyapunov spectra, plotted with a reduced 
index 1/2N on the abscissa; Lower panel: power spectra P i /2(h) for the modes corresponding to 
the smallest positive exponent, A; = 2at_3, of each spectrum. 

at the particle positions in the simulation cell for the 375-particle system with a density 
p = 0.4 familiar from Fig. El Naively, the upper surface would have to be classified as 
a transverse mode with the perturbation component 5y perpendicular to a wave vector k 
parallel to x, and the lower surface as a longitudinal mode with 5x parallel to k. Since the 
momentum perturbations are not nearly as parallel to the position perturbations (see Fig. 
EJ) as is (approximately) the case for hard disks, it is not sufficient to represent a mode by a 
single two-dimensional vector field such as (pq l \x,y) or tpp\x,y). Position and momentum 
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FIG. 11: Representation of a Lyapunov mode as periodic spatial patterns of the position pertur- 
bations 5xi and 5yi of the particles in the periodic box. The system consists of 375 WCA disks at 
a density p = 0.4 in a rectangular box with an aspect ratio A = 0.6. The Lyapunov spectrum for 
this system is given by the top-most curve of Fig. [21 The mode I = 746 is shown. 

perturbations need to be considered simultaneously for a characterization of the modes. The 
disappearance of the step structure in the Lyapunov spectrum and, hence, of the degeneracy 
for larger densities points to a complicated tangent-space dynamics we have not yet been 
able to unravel. We do not know of any theory which accounts for all of the numerical 
results presented in this paper. 
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